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Abstract 

We study the propagation of waves across fixed mesh refinement boundaries in linear 
and nonlinear model equations in 1-D and 2-D, and in the 3-D Einstein equations of 
general relativity. We demonstrate that using linear interpolation to set the data in 
guard cells leads to the production of reflected waves at the refinement boundaries. 
Implementing quadratic interpolation to fill the guard cells suppresses these spurious 
signals. 
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1 Introduction 



Wave propagation is an important phenomenon throughout all areas of physics, 
with applications typically involving multiple spatial and temporal scales. In 
numerical modeling of such problems, one strategy for dealing with the dis- 
parity in spatial and temporal scales is the use of a nonuniform or adaptive 
computational mesh. In this case waves can cross mesh refinement boundaries 
as they propagate through the computational domain. This paper focuses on 
interface conditions that will allow waves to travel smoothly across fixed re- 
finement boundaries, minimizing spurious reflections. 

The specific application that motivated this study is modeling the emission of 
gravitational waves from astrophysical sources such as binary black hole and 
neutron star coalescences. Such systems are among the most important sources 
for ground-based gravitational wave detectors such as LIGO and VIRGO [1,2], 
as well as the planned space-based LISA mission [3] . The gravitational waves 
produced typically have wavelengths ~ 10 — 100 times the scales of their 
sources Numerical simulations of such systems must therefore allow the signals 
to propagate from finely resolved regions around the sources into more coarsely 
resolved regions in the wave zones. Since the waveforms must be computed at 
large distances from their sources (i.e., effectivelly at infinity) for comparison 
with observations from gravitational wave detectors, the simulation domains 
must be made as large as possible. This can be achieved by incorporating 
several levels of successively coarser grids. 

The propagation of gravitational waves is governed by the Einstein equations, 
which are a coupled set of nonlinear partial differential equations [4]. These 
equations can be written in a variety of ways, but current practice in nu- 
merical relativity favors the use of the so-called BSSN formalism [5,6]. In this 
formalism, the Einstein equations are written as a system of quasi-linear equa- 
tions with first-order time derivatives and second-order spatial derivatives. In 
this paper we restrict our analysis to the "iterated Crank-Nicholson" update 
scheme, which is a second order accurate, explicit finite difference method that 
is currently in widespread use in the relativity community. It should also be 
noted that we consider mesh refinement only in space, not in time. In par- 
ticular, for our present analysis we use a common timestep across the entire 
computational domain. 

Adaptive mesh refinement (AMR) was first applied in numerical relativity to 
study critical phenomena in the 1-D collapse of a scalar field to form a black 
hole [7]. An early 3-D application focussed on evolving a single black hole 
[8]; this was followed by the use of fixed mesh refinement (FMR) to evolve 
a short part of a binary black hole evolution [9]. AMR was also employed to 
follow the propagation of gravitational waves through spacetime, first using 
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a single model equation that describes perturbations of a non-rotating black 
hole [10] and later in the 3-D Einstein equations [11], and to study inhomo- 
geneous cosmological models [12]. In these AMR studies the refinement and 
derefmement conditions were generally tuned so that the gravitational waves 
remained within the finely resolved regions. 

In this paper, we address the challenge of evolving gravitational wave signals 
across mesh refinement boundaries using FMR. Success in this endeavor is an 
essential component of gravitational wave source modeling, due to the large 
disparity in the scales of the sources and the waves. Our challenge amounts 
to choosing a prescription for coupling adjacent grid blocks when the blocks 
have different resolutions. Grid blocks are coupled through their guard cells, 
which must be filled using data from the blocks' interior cells. In hydrody- 
namics codes it is common practice to use a linear interpolation scheme for 
guard cell filling, with a possible adjustment for flux conservation across the 
interface between blocks [13,14,15]. We have found that this prescription is not 
adequate for the BSSN formulation of the Einstein equations. In particular, 
linear guard cell filling leads to unacceptably large reflections and distortions 
of the gravitational waves as they propagate from fine grid blocks to coarse grid 
blocks. Our solution to this problem is to use a guard cell filling procedure 
with quadratic-order accuracy orthogonal to the coarse-fine grid interface. 
The need for quadratic order guard cell filling has been previously demon- 
strated for elliptic boundary value problems with second order derivatives in 
[16,17]. With this prescription spurious wave reflections and distortions are 
reduced dramatically. 

Given the complexity of the full system of Einstein equations, we have chosen 
to analyze first a set of model wave equations in 1-D and 2-D that mimic some 
of the properties of the Einstein equations, as expressed in BSSN form.These 
simplified test beds have proved essential to understanding and correcting 
the problems that arise in the propagation of waves across mesh refinement 
boundaries. Since the solution we uncovered using these model equations has 
proved effective in curing the difficulties encountered in the Einstein equa- 
tions, we expect this work to be useful across a broad range of related wave 
propagation problems. 



2 Linear Wave Equation in 1-D: Evolution on a Uniform Grid 

The linear wave equation in 1-D is generally written in the form 

d 2 (f> _ d 2 (p 
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where <f> = <f>(x,t). Introducing the auxiliary variable H(x,t), we can cast 
Eq. (1) in a form that uses only first time derivatives: 



~dt 



d 2 



dx 2 



(2) 
(3) 



In this section we examine the system of equations (2)-(3) to understand 
the interface conditions needed for smooth propagation of waves across mesh 
refinement boundaries. In later sections, these conditions are applied to non- 
linear and multidimensional wave equations. 



2. 1 Discretization 



For the spatial discretization of equations (2)-(3), we take the data to be de- 
fined at the centers of the spatial grid cells and use standard O(Ax) 2 centered 
spatial differences [18]. To advance this system of ordinary differential equa- 
tions in time we use an O(At) 2 iterative method first suggested by M. Chop- 
tuik (see Ref. [19]). In the numerical relativity literature, this explicit update 
scheme is refered to as "iterated Crank-Nicholson". Each iteration has the 
form 

c +1 = # + At n, (4) 



where we use i to label the spatial grid, n to label the time steps, and <j>^ IT, 
to indicate intermediate values calculated during the iteration process. Note 
that the familiar Crank-Nicholson algorithm is obtained by setting <p i and Ilj 
equal to their time averages, + 0™)/2 and (II" +1 + II™)/2, respectively. 

For two iterations, the specific steps are as follows. Begin by applying the 
discretization (4)-(5) with <p = (fi n , II = IP to calculate a first approximation 
to <p n+1 and n n+1 : 

^ur 1 = €+^t it/ (6) 
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«n« +1 = n« + ^F(0«). (7) 



Average these new values with those at the starting time level n to get new 
values for and IT: 

(1) ^ = K (i vr i +^) (8) 

wn^K^n^ + n™). (9) 

Now perform a second iteration. Again applying (4)-(5) we find a second 
approximation to n+1 and n ra+1 : 

(2) 0™+! = 0" + At (^TTi (10) 

(2) n? +1 = n« + ^F(«0). (ii) 

Averaging again with the values at level n yields 

(2) ^ = |( (2) C +1 + C) (12) 
^n^K^n^ + n™). (13) 

A final update is carried out using these twice-iterated values: 

0«+! = ft + At ^TTi (14) 

nr +1 = nr + ^F(( 2 )0). (15) 

Clearly this algorithm can be carried out for any number of iterations. In the 
formal limit of an infinite number of iterations, it yields the usual Crank- 
Nicholson scheme. However, a von Neumann stability analysis shows that this 
iterative scheme is stable only when the number of iterations equals 2, 3, 6, 7, 
10, 11, etc, and the Courant condition At < Ax is satisfied. This was shown 
by Teukolsky [19] for the advection equation, but the conclusion holds as well 
for the wave equation in the form (2)-(3). Furthermore, the accuracy of the 
iterative scheme is second order for any number of iterations. We must carry 
out at least two iterations for stability, but continuing beyond two iterations 
does not reduce the truncation error. In this paper we follow the common 
current practice in numerical relativity and carry out precisely two iterations 
for our tests. 
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2.2 Evolutions on a Uniform Grid 

We first carried out uniform grid, or unigrid, evolutions of the discretized wave 
equation to provide a basis for comparison with mesh refinement runs. The 
initial data for <f> is taken to be a Gaussian wavepacket, 

<j>{x, t = 0) = A e~ x2/(T \ U(x, t = 0) = 0, (16) 

with A = 1 and a = 0.25. The spatial domain extends from x = — 4 to 
x = +4. Time evolution of this data produces two packets traveling with 
velocity v = ±1, each having amplitude A = 0.5 and the same value of o 
as the original packet. Here we will consider only the packet traveling to the 
right, in the region < x < 4. 

Figure 1 shows the evolution of this packet for two different resolutions. The 
coarser resolution is given by H — Ax = 0.045 (dotted line), which has ~ 10 
zones across the width of the packet at half its maximum amplitude. The 
solid line shows resolution h = H/2 = 0.0225. The time step is chosen to 
be At = Ax /A for a given spatial resolution, Ax. In the last few panels 
of Fig. 1 one can see a slight separation between the two curves. This is 
primarily due to numerical dispersion, which causes the phase velocities to 
deviate from unity. The phase velocity for a monochromatic wave propogating 
on a discrete, uniform grid is calculated in the Appendix, with the result 
displayed in Eq. (A. 11). According to this formula we expect the pulse (which 
has wavelength ~ 1) to propagate with speed ~ 0.999 on the fine grid and 
speed ~ 0.996 on the coarse grid. This translates into a separation between the 
two pulses of about 0.01 at time t = 3.37, which is the approximate separation 
seen in the last panel of Fig. 1. 

The time evolution of the absolute errors e = |0 ana i y ti C — 0numcricai| is shown 
in Fig. 2. The dotted line shows the errors en for the coarse resolution H, 
and the solid line is 4 x e h . Inspection of Fig. 2 shows that the two curves are 
nearly identical, demonstrating the second-order convergence of these runs. 
Note that the errors are approximately antisymmetric about the location of 
the pulse center. This is because the dominant source of numerical error is 
dispersion, which has the principle effect of shifting each wave pulse relative 
to the exact solution. 



3 Implementation of Mesh Refinement 

We use the Paramesh package [20] to implement the mesh refinement and 
parallelization in our codes. All of our codes use cell-centered data. Paramesh 



6 



works on logically Cartesian, or structured, grids and carries out mesh refine- 
ment on grid blocks. The underlying mesh refinement technique is similar to 
that of Ref. [21], in which grid blocks are bisected in each coordinate direction 
when refinement is needed. The grid blocks all have the same logical struc- 
ture, with nxb zones in the x— direction, and similarly for nyb and nzb. Thus, 
refinement of a block in 1-D yields two child blocks, each having nxb zones 
but with zone sizes a factor of two smaller than in the parent block. When 
needed, refinement can continue on the child blocks, with the restriction that 
the grid spacing can change only by a factor of two, or one refinement level, 
at any location in the spatial domain. Each grid block is surrounded by a 
number of guard cell layers that are used in computing finite difference spatial 
derivatives near the block's boundary. These guard cells must be filled using 
data from the interior cells of the given block and the adjacent block. 



Figure 3 shows a section of a 1-D grid in the vicinity of an interlevel boundary 
between two neighboring grid blocks. The fine grid covers the left half of 
the 1-D space, with cell-centered grid points labeled —1/2, —3/2, etc. The 
coarse grid covers the right half with cell-centered grid points labeled 1/2, 
3/2, etc. The fine and coarse blocks are offset from one another for clarity of 
presentation. One layer of guard cells is shown, with "G" marking the coarse 
grid guard cell and "g" the fine grid guard cell. These guard cells are filled 
with data from neighboring blocks or, if the block forms part of the edge of 
the computational domain, from appropriate outer boundary conditions. 



Paramesh can be used in applications requiring AMR, FMR, or a combina- 
tion of these. It handles the creation of grid blocks, and builds and maintains 
the data structures needed to track the spatial relationships between blocks. 
It takes care of all inter-block communications and keeps track of physical 
boundaries on which particular conditions are set, guaranteeing that the child 
blocks inherit this information from the parent blocks. In a parallel envi- 
ronment, Paramesh distributes the blocks among the available processors to 
achieve load balance, maximize block locality, and minimize inter-processor 
communications. 



For the work described in this paper, we are using FMR. For simplicity, we 
use the same timestep, chosen for stability on the finest grid, over the entire 
computational domain. At the mesh refinement boundaries, we use a single 
layer of guard cells as shown in Fig. 3; special attention is paid to the restriction 
(transfer of data from fine to coarse grids) and prolongation (coarse to fine) 
operations used to set the data in these guard cells, as discussed in the next 
subsection. 
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4 Linear Wave Equation in 1-D: Evolutions with Fixed Mesh Re- 
finement 

We now carry out evolutions of 1-D linear waves that encounter a change in the 
grid resolution at a fixed location. For the gravitational wave applications in 
which we are interested, waves will be generated in a finely resolved region and 
then travel out into more coarsely resolved regions. We thus start our initial 
wave packet, given by Eq. (16), in a region of fine resolution h = 0.0225 around 
the origin. The spatial domain is again — 4 < x < 4. As before, the initial 
wave packet splits into two identical packets traveling in opposite directions. 
Each of these packets then encounters a fixed refinement boundary, located 
at x = ±2.1, and crosses into a region of coarser resolution H = 2h. In the 
following discussions, we focus only on the region x > 0. 

We first use the default Paramesh linear interpolation to set the value of the 
data in the guard cells on both the coarse and fine grids. With this prescription 
for guard cell filling, the coarse grid guard cell value of any function / is given 
by linear interpolation, 

/G = ^(/-S/2 + /-l/2). (17) 

The value of / in the fine grid guard cell "g" is then given by a linear interpo- 
lation using coarse grid values, / g = (/ G + 3/i/2)/4. Combined with Eq. (17), 
this gives 

/g=^(/-S/2 + /-l/2 + 6/1/2). (18) 

Note that this guard cell filling (GCF) procedure uses the points / G and /i/ 2 
on the coarse grid to obtain / g ; this is in contrast to the direct approach, which 
uses the nearest points f-1/2 and /1/2 (cf. Eq. (23)). The prescription (17)- 
(18) for GCF has errors of order h 2 and is the default linear GCF method 
in Paramesh. We will refer to this procedure as linear GCF in this paper. 
The results of using linear GCF are displayed in Fig. 4, which shows the time 
evolution of the absolute errors e. The dotted line shows the run with linear in- 
terpolation at the interface boundary, and the solid line the results of a unigrid 
run at the fine grid resolution. As the packet passes through this boundary, 
a reflected wave is generated propagating to the left. The transmitted wave 
continues traveling to the right into the coarse grid region. 

In large scale simulations of the Einstein equations with several levels of re- 
finement, such spurious reflected waves can seriously degrade the quality of 
the results. Globally increasing the resolution until the reflected waves reach 
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acceptably small amplitudes is generally not possible in 3-D. We thus need a 
better way to control the behavior of the signals crossing the interfaces. 

To this end, we implemented direct quadratic interpolation (i.e., using the 
nearest 3 data points) to set the data in the coarse and fine grid guard cells. 
Refer again to Fig. 3. For the fine grid guardcell "g" , quadratic interpolation 
yields [18] 

/ g = ^(-3/- 3 /2 + 10/_i /2 + 8/i /2 ). (19) 



The coarse grid guard cell "G" is filled by matching first derivatives across 
the interface, 

fl/2 - /g _ 7g - f-1/2 

H ~ h ' 1 Uj 



where H = 2h. This step, which ensures that the solution is smooth across 
the interface, can be viewed as "flux matching" where the gradient of / plays 
the role of the flux. By combining the derivative matching condition with the 
formula for f g we find 

/g = ^(6/_3/2 + 10/_i /2 -/i /2 ). (21) 



This same result for /g can be obtained by direct quadratic interpolation. 
These formulae for GCF have errors of order h 3 . 

The absolute errors obtained when using quadratic interpolation are shown as 
the dashed line in Fig. 4. Note that the reflected wave has been greatly reduced. 
Additional simulations, in which the size of the zones is everywhere decreased 
by successive factors of two, show that with quadratic GCF the code is second- 
order convergent. On the other hand, with linear GCF, the reflected pulse is 
first-order convergent. The transmitted pulse also aquires first-order errors at 
the interface with linear GCF. As the transmitted wave propagates through 
the coarse grid region, second-order errors due to dispersion and dissipation 
eventually dominate over the first-order errors introduced at the interface. At 
that point, the transmitted pulse can appear second-order convergent. 

We also conducted tests using a one-dimensional periodic domain consisting of 
20% fine grid and 80% coarse grid. A wave pulse was allowed to cycle through 
the domain multiple times. These tests clearly show that with quadratic guard 
cell filling, but not with linear guard cell filling, the code is second-order 
convergent. We also used this test code to check the stability of the interface 
conditions. After thousands of cycles of the wave pulse through the refined 
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region, there were no signs of instability with either linear or quadratic guard 
cell filling. 

In the appendix we present a detailed analytic treatment of wave propagation 
across mesh refinement boundaries that complements our numerical experi- 
ments. There we compute the reflection coefficient R and transmission coeffi- 
cient T for a monochromatic (single frequency) wave traveling on a grid with 
fixed mesh refinement, for various methods of GCF. The wave travels from a 
fine grid region with resolution h into a coarse grid region with resolution 2h. 
Figure 5 shows the absolute value of the reflection coefficient |R| for linear 
GCF (17)-(18) (dashed curve) and quadratic GCF (21)-(19) (solid curve). 
The dotted curve shows the results for direct linear interpolation, defined by 

/o = i(/-s/2 + /-i/2). (22) 

for the coarse grid guard cell and 

/g = ^(/-i/2 + 2/i /2 ). (23) 



for the fine grid guard cell. Direct linear interpolation, like the default linear 
GCF in Paramesh, has errors of order h 2 . The curves of Fig. 5 are plotted as 
functions of the wavelength in the fine grid region divided by the fine grid cell 
size h. Equivalently, we can interpret the horizontal-axis values as the number 
of fine grid points per wavelength. 

For our 1-D wave equation tests, the Gaussian packet behaves roughly like 
a wave of wavelength A ~ 1. With h = 0.0225, this corresponds to about 
X/h = 44 fine grid points per wavelength. From Fig. 5 we see that the reflection 
coefficient for linear interpolation is about |R| = 0.02 while that for quadratic 
GCF is |R| = 0.0003. With an incident pulse amplitude of 0.5, we expect a 
reflected wave amplitude of about 0.01 for linear GCF and less than 0.0002 
for quadratic GCF. This reflected pulse for the linear case is clearly seen in 
Fig. 4. 

The importance of minimizing spurious reflections from grid interfaces has 
been emphasized above. It is equally important to minimize the distortion of 
waves that pass through a grid interface. The errors in the transmitted wave 
pulse for linear and quadratic GCF are shown in the region x > 2.1 of the last 
few panels of Fig. 4. Note that the errors for quadratic GCF are actually larger 
than the errors for linear GCF. This surprising result is explained as follows. 
Observe that the errors for the two fixed mesh refinement simulations, as well 
as for the unigrid run (solid curve), are approximately antisymmetric about 
the pulse center. The errors in each case, as in the unigrid tests discussed in 
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Section 2, are primarily due to dispersion. Dispersion causes the wave pulses 
to fall behind the exact solution during propagation, giving rise to the errors 
shown in Fig. 4. This effect is greater for the two runs with fixed mesh re- 
finement because, beyond x = 2.1, the grid resolution is lower than for the 
unigrid run. However, with mesh refinement, the transmitted pulse will also 
suffer a phase error which has the effect of artificially shifting the pulse along 
the x-axis. In the case of linear GCF, there is a relatively large positive phase 
error in the transmitted wave. This phase shift partially compensates for the 
negative shift caused by dispersion. As a result the size of the largest peaks 
in the error for the transmitted wave, for the particular test shown in Fig. 4, 
is smaller with linear GCF than with quadratic GCF. 

Figures 6 and 7 show the absolute value of the transmission coefficient |T| 
and the phase of the transmission coefficient ip = arctan($5(T)/K(T)) for a 
monochromatic wave, for linear, quadratic, and direct linear interpolation. 
These graphs are obtained from the analysis in the Appendix. From Fig. 6 
it is clear that at any wavelength (any resolution) the error in amplitude for 
the transmitted wave is smaller for quadratic GCF than for linear GCF. 2 
The dominant source of error for the transmitted wave is actually phase error, 
shown in Fig. 7. The magnitude of this error for quadratic GCF is much 
smaller than that for linear GCF. For a wavelength of A ~ 1, the linear guard 
cell filling produces a phase shift of about ip = 0.024, while quadratic GCF 
gives a phase shift of about p = —0.00028. For the tests shown in Fig. 4, the 
positive phase for linear GCF translates into a shift along the positive 
of about Sx = \<p/(2ir) ~ 0.004. With quadratic GCF, the pulse is shifted in 
the negative direction, but by a much smaller amount Sx ~ —0.00004. Close 
inspection of the data for the two transmitted pulses shows that they indeed 
have a separation of 5x m 0.004. For linear GCF, this phase shift pushes 
the wave pulse forward and artificially compensates for the phase lag caused 
by dispersion. In general, there is no reason to expect the cumulative phase 
lag due to dispersion to be close in magnitude (but opposite in sign) to the 
phase advance caused by transmission through various grid interfaces. Thus, 
the relatively small transmission error seen in Fig. 4 for linear GCF should be 
viewed as an accident of the particular example, not a generic result. 



2 At low resolution, that is, for wavelengths less than about 28h, direct linear GCF 
has the smallest error for the transmitted wave amplitude. However, as discussed 
in the Appendix, as the resolution is increased |T| is much closer to 1 for quadratic 
GCF. Also note from Fig. 7 that direct linear GCF has large phase errors for the 
transmitted wave. 
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5 Nonlinear Wave Equation in 1-D 



The next step in developing model equations to test these interface conditions 
is to add nonlinear terms similar to those found in the Einstein equations. 
This produces the following nonlinear wave equation 



52J, /o,\2 /o,\2 



where d and e are arbitrary contants. Again introducing the auxiliary variable 
n(x, i), we get the first order system 

£-n (25) 



Using the discretization introduced in § 2.1, we have 

= € + (At)n, (27) 



+ e(At)(^il) 2 . (28) 

Equations (27) and (28) are updated following the steps given in (6)-(15). 

We consider the case d — — e — 1 and set up an initial Gaussian wave packet 
centered on the origin using the prescription given by Eq. (16), with H(x,t = 
0) = 0. This splits into two identical packets traveling in opposite directions, 
each having amplitude A = 0.38 and width a = 0.25. We use the spatial 
domain —4 < x < 4 and set fixed refinement boundaries at x — ±2.1. The 
fine grid around the origin has resolution h = 0.0225 and the coarse grid 
regions have resolution H = 2h. We focus on the region x > 0. 

The results are shown in Figure 8. Since we do not have an analytic solution for 
Eq. (24), we display the actual solution and use unigrid runs for comparison. 
In addition, the vertical scale is chosen to zoom in on the region around the 
base of the packet (i.e., near = 0), where the differences between the runs are 
the most apparent. The thin solid line shows the solution for a unigrid run at 
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the coarse resolution H, and the thick solid line shows a unigrid run at the fine 
resolution h. Runs in which the packet encounters a refinement boundary are 
shown using a dotted line (linear GCF) and a dashed line (quadratic GCF). 
As we saw before, a reflected wave is generated when the packet crosses the 
refinement boundary using linear GCF; these effects are much less noticeable 
when using quadratic GCF. As in the case of the linear wave equation, the 
code is second-order convergent when using quadratic GCF. 



6 Wave Equation in 2-D 



As a next step, we consider the wave equation in 2-D. The evolution of cylin- 
drically symmetric waves on a 2-D Cartesian mesh provides an ideal test 
problem in which the signals cross mesh refinement boundaries that are, in 
general, not perpendicular to their directions of propagation. 

The 2-D model wave equation takes the form 

9 2 d 2 <p d 2 <p (d<p\ 2 (d<P\\ (d<p\ 2 



where d, e±, ei are contants. With the auxiliary variable H(x, y, t), we can write 
this in a form using only first-order time derivatives: 

£ = n (30, 
du d 2 <p d 2 <p 2 ( d <t>\\ (WY 



Using the discretization introduced in § 2.1, we have 

<Pt l = + ( A *)% (32) 



+ rf(At)(n,,) 2 + e 1 (At) (^'• : -^ ; 1 
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As before, Eqs. (32) and (33) are updated following the steps given in Eqs. (6)- 
(15). 

In this section we consider two types of GCF, the default Paramesh linear order 
GCF and a quadratic GCF scheme. The linear GCF is depicted in Fig. 9. First, 
each coarse grid guard cell (open diamond) is filled as a linear combination of 
the surrounding fine grid points (solid circles). The fine grid guard cells (open 
circles) are then filled using a linear combination of the surrounding coarse 
grid points (open and solid diamonds). 

Quadratic GCF is depicted in Fig. 10. As a first step, the coarse grid guard 
cells are filled from a linear combination of the four surrounding fine grid 
guard cells. These values are only used at fine grid corners, and will soon 
be overwritten. Linear interpolation of the coarse grid cells (solid diamonds) 
parallel to the coarse-fine interface is used to compute intermediate values 
marked with open boxes in Fig. 10. These intermediate values, along with the 
two fine grid cells (solid circles) directly across the interface, are then used to 
obtain a quadratic fit for the fine grid guard cells marked with open circles. 

Finally, as in the 1-D case, the coarse grid guard cells are filled by "flux 
matching", that is, matching derivatives across the interface. Specifically, we 
consider the first derivative at the midpoint of Fig. 10, that is, at the point 
midway between the coarse grid guard cell (open diamond) and the interior cell 
directly across the interface (closed diamond). Derivative matching consists in 
equating the first derivative computed from these coarse grid cells with the 
second-order accurate first derivative obtained from the four fine grid cells 
(open and closed circles) that surround the midpoint. 

The algorithm described here for quadratic GCF is similar to the one de- 
scribed by Martin and Cartwright [23]. The main difference is that we use 
linear interpolation of coarse grid values parallel to the interface to obtain in- 
termediate values (the open boxes in Fig. 10), whereas Martin and Cartwright 
use quadratic interpolation. Also note that our algorithm can be applied with- 
out modification at fine grid corners, where the corner of a coarse grid block 
is surrounded by fine grid blocks. Recall that in the first step, coarse grid 
guard cells are filled by linear restriction from the fine grid. This allows the 
interpolation parallel to the interfaces to be carried out without the use of 
one-sided extrapolation. Finally, we point out that GCF at fine grid corners 
is ambiguous, since there are different ways to deal with them; either of the 
two coarse-fine interfaces that intersect at the corner can be used or intepola- 
tion using a stencil diagonal to the interfaces can also be used. Note that only 
mixed derivatives are affected by the corners when using centered differencing. 
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In our code we do not treat the corners as special. At a corner our code nat- 
urally selects one of the two interfaces and carries out a linear interpolation 
parallel to that face to obtain intermediate values. 

The initial data for our tests is taken to be a cylindrically symmetric wavepacket 
centered on the origin, with 

<f>(x,y,t = 0) = Ae^ x2+y2)/a2 (34) 



and n(x, y, t — 0) = 0. We choose the amplitude A — 1, and the width of the 
pulse by a = 0.25. Quadrant symmetry is imposed by using mirror-symmetry 
boundary conditions along x = and y = 0. The computational domain then 
covers the region < x < 4.3125 and < y < 4.3125. This packet is initially 
confined to a fine grid region of resolution h = 0.0225. As the packet expands, 
the wavefront crosses a fixed mesh refinement boundary into a region of coarser 
resolution H = 2h. 

Setting = ei = e2 = in Eqs. (30) and (31) allows this packet to evolve 
under a linear equation. Figure 11 shows the results of using quadratic GCF 
to set the values of the guard cell data. Here, is shown at four consecutive 
times. The expanding wavefront encounters mesh refinement boundaries at 
x — 2.1 along the a;— axis and at y — 2.1 along the y— axis. Note that the wave 
passes smoothly across the interface. 

A comparison of unigrid and fixed refinement runs is shown in Fig. 12. Here, 
is plotted along a portion of the x— axis at a fixed time. The unigrid run 
(solid line) at the fine grid resolution shows the extended "tail" of the out- 
going cylindrical wave front. The run with linear GCF (dotted line) shows 
a reflected wave traveling back into the fine grid region as the wave passes 
through the refinement boundary. In the run with quadratic GCF (dashed 
line), this spurious signal has been nearly eliminated. 

Similar results are achieved when this wave packet is evolved according to a 
nonlinear equation, d = —e\ = — e 2 = 1. The structure of the grid and loca- 
tion of the refinement boundary are the same as for the 2-D linear equation. 
Figure 13 displays the results of along the at a fixed time. Notice 

that the run with linear GCF (dotted line) shows a significant reflected wave. 
In contrast, the run with quadratic GCF (dashed line) is close to the one with 
a uniform grid (solid line). 



15 



7 The Einstein Equations in 3-D 



We are now ready to apply the techniques developed in our model equations 
to the propagation of gravitational waves in 3-D, which is governed by the 
vacuum (or source-free) Einstein equations. We write these equations in terms 
of the "3 + 1" spacetime split [4], in which the initial data is specified on some 
3-D spacelike slice and then evolved forward in time. Within this framework, 
the metric takes the form 

ds 2 = -a 2 dt 2 + gij (dx i + (?dt) (dx j + (3 j dt) . (35) 



We use units in which both the speed of light c = 1 and the gravitational 
constant G — 1. Lowercase Latin letters are used to denote spatial indices, so 
that i, j = 1, 2, 3. To simplify the notation throughout this section, we use the 
summation convention: if any expression has one index as a superscript and 
the same index as a subscript, summation over all values that index can take 
is implied [22]. The geometry of the given spacelike slice is described by the 
3-metric g^. The lapse function a governs the advance of proper time across 
the surface, and the shift vector f3 % the motion of the spatial coordinates within 
the hypersurface as the data is evolved forward in time. Both a and f3 l are 
freely-specifiable functions of space and time; for the rest of this section, we 
use the choice a — 1 and f3 l = 0. 

In the standard ADM spacetime split [4], the Einstein equations can be written 
in terms of and the extrinsic curvature of the hypersurface where 



Following current practice in numerical relativity, we use the BSSN formalism 
[5,6] in which the Einstein equations are written in terms of conformal variables 
{tp, K, g^, Aij, P} defined as follows: 



e 4 ^ = det(^) 1 / 3 (37) 

~ 9ij = e~^g l3 (38) 

K ee gVKn (39) 

4- ee e -*(tfy - l - 9ij K) (40) 

r ee -djF . (41) 
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Here g l i is the inverse of the conformal metric gij. We use the notation dj = 
d/dxi for spatial derivatives. 

In terms of these conformal variables, with the gauge choices a = 1 and f3 l = 0, 
the vacuum Einstein equations become 

dib 1 

it = s aK < 42 > 



dt 

dK 1 



-24- (43) 
+ AijA? j (44) 



9* 3 

^ = i?J F + - 2A,^ (45) 
r?f* 2 

— = 2(P jfc ^ - + 6^9^). (46) 



Here, T*^ are the connection coefficients associated with g^, defined by 

f k ji = ^g mk (d z g mj + dfg mi - d m ~g 3l ) , (47) 



and 

~ A H = f~ g j*A lk , A l j = f~A ty (48) 

The superscript "TF" denotes the trace-free part of a tensor, so that Rj? = 
Rij — gijR/3, where R = g mk R m k- The Ricci curvature tensor Rij is defined 
by 

p fj "pfc _ f) -pk _|_ pfc pm _ pfe pm //|Q\ 

Although the set of equations (42)-(46) is considerably more complicated than 
our model equations, there are notable similarities. In particular, the conformal 
metric g^ plays the role of the function 0, while Aij takes the role of n. Looking 
at (47) and (49), we also see that the Rj^ term in Eq. (45) contains second 
spatial derivatives of g^. 

The lessons learned from the model equations in 1-D and 2-D can be applied 
successfully to the Einstein equations in 3-D, as we demonstrate by evolving a 
weak gravitational wave. We use the analytic solution to the linearized Einstein 
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equations found by Teukolsky [24]; since this is given in closed form, we can 
then compare the numerical results directly with this analytic solution. We 
choose the even parity, L = 2, M = solution, which is given by 



ds 2 = -dt 2 + (1 + Af rr )dr 2 + (2Bf r . e )rdrd9 

+(2Bf r4> )r sin Odrd<j> + (1 + Cf$ + Af$)r 2 d9 2 
+[2(A- 2C)f 6<p }r 2 sin 9d9d<f> 

+ (1 + C •f$}+Af$y sin 2 9d<j> 2 . (50) 

Here, 



A = 3 
B = - 



—it- + t- + —r 



F(3) 3F( 2 ) 6F« 6F 



i?(4) 9F( 2 ) 21F^ 21F 



+ 



+ 



F = F(t-r), F (n) = 



+ 



+ 



dx T 



(51) 
(52) 
(53) 
(54) 



- x=t—r 

where F is a generating function. We use the form 



F(x) = ^e"^ 2 , (55) 
or 

with two free parameters, A and u. Here we have specified an outgoing wave 
solution F = F(t — r); an ingoing wave solution can be obtained by using 
F = F(t + r). 

For this even-parity, M = case, the angular functions f\j are: 



f rr = 2-3 sin 2 9 (56) 

f r9 = -3sm9cos9 (57) 

U = (58) 

/ff=3sin 2 ^ (59) 

A 2) = -l (60) 

fe^ = (61) 

/# = "iff (62) 

/# = 3 sin 2 9 - 1. (63) 
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We present results for a gravitational wave crossing two fixed mesh refinement 
boundaries into regions with successively coarser resolution. We start with a 
wave packet composed of a linear combination of one initially ingoing and one 
outgoing wave, each having amplitude A = 1CT 6 and width uj — 1. This packet 
is centered on the origin in a fine grid region of resolution h = 0.0416667. 
The successively coarser regions have resolutions 2h and Ah, with the first 
refinement boundary at r = y/x 2 + y 2 + z 2 = 4.5 and the second at r = 9.0. 
To complete the initial data we take = so that K — and = 0. 
Octant symmetry is imposed by using mirror-symmetry boundary conditions 
along x = 0, y = 0, and z = 0. The computational domain covers the regions 
< x < 12 and similarly for y and z. 

As the evolution proceeds, the outgoing waves travel directly toward the outer 
boundary of the grid. The initially ingoing waves first travel toward the origin, 
then reflect and move outward. As the overall signal propagates outward, it 
leaves flat spacetime behind. 

Figure 14 shows the evolution of these waves when linear GCF is used. The 
function g zz — 1 is plotted as a function of x and y in the z = plane at 4 
successive times. Note the presence of spurious reflected signals as the waves 
pass through the fixed mesh boundaries. These problems are greatly reduced 
when quadratic GCF is used, as shown in Fig. 15. A comparison of runs with 
linear (dotted line) and quadratic (dashed line) GCF and the analytic solution 
(solid line) is shown in Fig. 16. The reflected waves are essentially eliminated 
by the use of quadratic GCF. 

Finally, Fig. 17 demonstrates the second-order convergence of the code by 
comparing the results of the run in Fig. 15 with a run that differs only by 
having the size of the grid zones a factor of 2 larger throughout. Both runs 
use quadratic GCF. The L2 norm of the absolute error e is calculated over 
each simulation domain, and plotted as a function of time. The solid triangles 
connected by the solid line show e for the run in Fig. 15, and the filled boxes 
connected by the dotted line show the errors for the lower resolution run 
multiplied by 4. 



8 Summary 

We have examined the propagation of waves across fixed mesh refinement 
boundaries, starting with simplified linear and nonlinear model equations in 
1-D and 2-D, and progressing to the 3-D Einstein equations of general rela- 
tivity. The numerical evolutions were carried out using centered spatial differ- 
ences and the explicit iterated Crank-Nicholson time update method, giving 
second-order accuracy. Our results show that using linear GCF produces spu- 
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rious reflected waves as the signals cross refinement boundaries, and that these 
are greatly suppressed by using quadratic GCF. In particular, quadratic GCF 
preserves the second-order convergence of the numerical evolutions. Our nu- 
merical results are complemented by a detailed analytic treatment of waves 
crosing refinement boundaries in 1-D in the Appendix. 

While quadratic GCF is straightforward to describe and implement in 1-D, the 
situation becomes more complicated in 2-D. In particular, intermediate values 
parallel to the mesh refinement interface must be calculated in the 2-D case. 
We have found that using linear interpolation to obtain these intermediate 
values, combined with quadratic interpolation for the final values, maintains 
the second-order convergence. The procedure used for quadratic GCF in 2-D 
generalizes in a straightforward manner to the 3-D case. 

The techniques presented here appear to be robust in the sense that they con- 
tinue to produce excellent results as our test problems increase in complexity. 
Quadratic GCF successfully eliminates most of the spurious reflected waves in 
both linear and nonlinear model equations in 1-D and 2-D. The 3-D Einstein 
equations present a much larger and more complex system of equations. In the 
test case presented here, the evolution of a weak gravitational wave, quadratic 
GCF continues to perform well, even as the signals cross two successive mesh 
refinement boundaries. We fully expect that these techniques will also yield 
excellent results for strong gravitational waves, which activate the nonlinear 
terms in the Einstein equations. Such evolutions require various technical dif- 
ferences in the gauge choices (a and f3 l ) as well as in the formulation of the 
initial data. We are currently working on such models, and will report on them 
elsewhere. 
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A Appendix: Analysis of Numerical Wave Propagation in 1-D 

In this appendix, we present a more detailed analysis of the propagation of 
linear waves in 1-D with the discretization described in Sec. 2. We begin by 
deriving some basic results for uniform grids (Sec. A.l) and follow with a study 
of wave propagation across a fixed mesh refinement boundary (Sec. A. 2). Here, 
we do not address the issue of instabilities that might arise due to coupling 
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between fine and coarse meshes [25]. However, as noted in Sec. 4, our numerical 
tests show no signs of instability. 



A.l Wave Propagation on a Uniform Mesh 



We will employ matrix notation to facilitate the analysis in this Appendix. 
First, we collect the field variables 0, II into the column vector 



f=(5)- (A-l) 



Equations (2)-(3) can now be written as 

~dt = (<9 2 /L 2 l) V ■ (A ' 2) 



As usual VJ 1 will denote the vector of grid functions at timestep n and grid 
point j. 

The iterated Crank-Nicholson method described in Sec. 2.1 is built from suc- 
cessive applications of the basic operator 

I) , (A.3) 



where d 2 VJ l = {Vj l +1 — 2VJ 1 + V™_ 1 )/Ax 2 . With two iterations, the update of 
the variables V n by one full timestep is accomplished by the operator 



M = I + AtQ 



2 ^ V 2 



(A.4) 



The stability, dissipation, and dispersion properties are obtained by consider- 
ing discrete plane wave solutions, 

V n = We iu}nAt e~ ikjAx , (A.5) 



where W is a constant vector (independent of n and j). Inserting this ansatz 
into the update equation V^ n+1 = MVJ 1 , we find 

e**"W-( 1_2A2 ^-^W fA6) 

6 W ~ V -4A 2 (1-A 2 )/Ai 1-2A 2 ) W ' [A - b) 
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where 



A = ^ sin(A;A:r/2) . (A.7) 



Thus, W is an eigenvector with eigenvalue e ' for the matrix that appears 
in Eq. (A. 6). The eigenvalues are obtained in the usual way with the result 
e iuAt _ y — 2 A 2 ± 2«A(1 — A 2 ). This is the dispersion relation giving the 
complex frequency w as a function of wave number k. We can, without loss 
of generality, consider only plane wave solutions (A. 5) with positive frequency 
£ > 0, where £ = K(cj) is the real part of oo. Then the ± sign in the dispersion 
relation must be set equal to the sign of the wave number k. The dispersion 
relation then becomes 

e ^At = i _ 2A 2 + 2i|A|(l - A 2 ) (A.8) 



and £ = dt(u) is positive. The eigenvectors W corresponding to these eigen- 
values are straightforward to compute. Choosing the first component of W to 
be unity, we find 



We note for later reference that We tk i Ax is an eigenvector of the basic oper- 
ator Q with eigenvalue 2i\A\/At. 

The finite difference scheme is unstable if the magnitude of the amplification 
factor, \e iujAt \, is greater than unity. From Eq. (A.8) we find that \ e iujAt \ 2 < 1 
implies A 2 < 1. This inequality will be satisfied for all wave numbers k only if 
At < Ax. This is the Courant limitation on the timestep for the wave equation 
(2)-(3) discretized with twice-iterated Crank-Nicholson. 

The phase velocity is found from the real part of the frequency £ = $t(u). 
From the dispersion relation (A.8), we find 

:..M = arcsin I _W^£^ = | ( A. 1 ) 



yi -4A 4 (1 - A 2 ) 



The phase velocity is then 

C(A )4 = ^LL (a. ii) 

k 2-na Ax K J 
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where a = At/ Ax is the Courant factor and A = 27r/k is the wavelength. The 
dissipation is found from the imaginary part of the frequency, rj = Of (a;). Since 
the wave amplitude varies like ~ e~ r,nAt , we see that the amplitude drops 
by a factor 

e - V At = | e iu,At| = ^/i _ 4A 4 (1 - A 2 ) . (A.12) 



for each timestep. 



A. 2 Wave Propagation with FMR 



Now consider a two-level refined mesh, with fine grid Axf on the left and 
coarse grid Ax c on the right. We will assume that the refinement jumps by a 
factor of 2, that is, Ax c = 2Axf. The mesh will be labeled as shown in Fig. 3. 
Thus, VT1/2? V-3/2 > e ^ c - are the fine grid functions and Vy 2 , ^3/2 > e ^ c - are the 
coarse grid functions. 

As a first step towards analyzing the wave reflection and transmission at the 
interface, we relate the wave numbers in the coarse and fine grid regions. 
Consider a monochromatic solution that varies like <f> ~ e t£,nAt across the 
entire mesh. Specifically, we assume that the coarse and fine grid frequencies 
are the same, £ c = £/, and that the coarse and fine grid time steps are the 
same, At c = Atf. From the dispersion relation, Eq. (A. 8), we can compute 
tan(£At) in the coarse and fine grid regions and equate the results: 

2|A C |(1 - Ag) _ 2|A / |(l-Af) 

1-2A2 " l-2Af • [AA6) 



Here, A c = (At/Ax c ) sm(k c Ax c /2) and similarly for Aj. This relation has the 
form /(|A C |) = /(|A/|) where |A C | and |A/| vary between and 1. It is easy to 
show that the function /(|A|) is monotonic and therefore invertible. It follows 
that the only solution of Eq. (A. 13) is 

|A C | = |A/| . (A.14) 



This equation shows that the coarse and fine grid wave numbers k c and kf are 
related by 



K = ± 



Ax, 



arcsm 



Ax c 
Ax f 



sm(kfAxf/2) 



(A.15) 
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The two cases corresponding to the ± sign indicate that the wave propagation 
direction on the coarse and fine sides of the interface need not match. Thus, 
we can have a right moving wave in the coarse grid region connected to both 
right moving and left moving waves in the fine grid region. 

From the result (A. 14) we see that the rate of dissipation (A. 12) of a wave, 
governed by r] = $>(u), is the same in the coarse and fine grid regions. We also 
see that the relative phase between the two components and II of the wave, 
Eq. (A. 9), is the same in coarse and fine regions. The phase velocity (A. 11), 
and hence the amount of dispersion, differ in the coarse and fine grid regions, 
since the wave numbers k c and kf are not equal. 

According to Eq. (A. 15) k c is real only for kfAxf < tt/3, that is, for Xf/Axf > 
6. If \f/Axf < 6, then k c is complex and the plane wave solution (A. 5) will 
contain a spatial dependence in the coarse grid region that is either expo- 
nentially damped or grows exponentially. Note that, although k c might be 
complex, A c is real (assuming kf is real) and equal to ±Aj. It follows that, 
whether k c is real or complex, the Courant stability condition | e ^ A *| 2 < 1 is 
satisfied in the coarse grid region if it is satisfied in the fine grid region. 

For the remainder of this appendix we will focus on the case of practical 
interest, where Xf/Axf > 6 and k c is real. The plots in Figures 5, 6, and 7 
have been restricted to Xf/Axf > 10 for clarity of presentation. Each of the 
curves in those plots reaches a finite value at Xf/Axf = 6. 



A. 2.1 Matching solutions 

At this point we have shown that waves of frequency £ have wave number k c 
on a coarse grid, wave number kf on a fine grid, and that these values are 
related as in Eq. (A. 15). We will now construct a solution with frequency £ 
that spans the entire non-uniform grid. To begin, consider the vector 

C W ( e -ik f jAx f + Re ik f jAxA ■ < q 

Vi = { > . . \ ' (A.16) 

3 [ W (Te-* k ^ Ax c) , j > . 



We will show that for an appropriate choice of the coefficients R and T the 
vector Vj is an eigenvector of the basic operator Q with eigenvalue 2i\A\/At. 
For points away from the interface, namely, the points j < —3/2 and j > 
3/2, this conclusion follows from our earlier observation that on a uniform 
grid We ±ikjAx is an eigenvector of Q with eigenvalue 2i\A\/At. The same 
argument cannot be applied to the points 1/2 and —1/2 surrounding the 
interface because the stencil for the discrete derivative operator d 2 appearing 
in Q extends across the interface. Thus, when computing d 2 Vj for j = ±1/2, 
we must use guard cell information. 
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In the main text we discussed various choices for guard cell filling, such as the 
Paramesh linear GCF of Eqs. (17)-(18) and the quadratic GCF of Eqs. (19)- 
(21). For the purpose of presenting the analysis, we will focus instead on the 
direct linear GCF of Eqs. (22)-(23). In the present notation, these relations 
are 

VS = ^- 3/2 + V? 1/2 ), V g n = i(^ 1/2 + 2Vfa) (A.17) 



Now, for grid points that are not adjacent to the interface, the operator d 2 
takes the usual form, 

° V i I (V- - 2VJ 1 + V")/Axl , j > 3/2 . [AA * } 



But for grid points adjacent to the interface, d 2 must use guard cell values 
given by (A.17). Consequently, we find 



d 2 V\ /2 = (V g n - 2V? 1/2 + V^ /2 )/Ax) 

= (2V l n /2 - 5VT 1/2 + 3VT 3/2 )/(3A4) (A.19) 

and 



&K,2 = (VST/2 - 2vr /2 + v£)/Ax 2 c 

= (2V;) 2 - W? /2 + Vl\ /2 + V^ /2 )/(2Ax 2 c ) . (A.20) 
for d 2 acting at grid points j = ±1/2. 

We now impose the requirement that the vector Vj of Eq. (A. 16) is an eigenvec- 
tor of Q with eigenvalue 2i\A\/At at the points j = ±1/2 adjacent to the in- 
terface. Using the discretization (A.19) the relation QV-i/ 2 = (2i\A\/ At)V-±/ 2 
yields 

i(20 1/2 - 50_ 1/2 + 30_ 3/2 ) = 0-i/2 • (A.21) 



Here, <pj is the first component of the ansatz vector Vj. Similarly, with the 
discrete operator (A.20), we find that QV\/ 2 = (2i\A\/ At)Vi/ 2 implies 



. , , , /2i|AI N 



■ (203 /2 - 40i /2 + 0-1/2 + 0-3/2) = 01/2 • (A.22) 
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These two equations can be solved for the two coefficients R and T. The result 
is 



_ 3E 2 E} - Ej - ElE) - Ej 
l + E} + E 2 E 2 f -3E 2 Ej ' 

E c (3 + 2E}-2Ej-3Ej)/(2E f ) 

1 + Ej + E 2 E 2 — 3E 2 Ej ' 1 ' 

where the shorthand notation 

E c = e ikcAxc/2 , E f = e ik f Ax f /2 (A.24) 



has been used. 



At this point we have succeeded in showing that the vector Vj of Eq. (A. 16), 
with R and T chosen as in Eqs. (A. 23), is an eigenvector of Q on the non- 
uniform grid. The corresponding eigenvalue is 2i\A\/At. A short calculation 
shows that Vj is also an eigenvector for M, Eq. (A. 4), with eigenvalue 1 — 
2A 2 + 2i|A|(l — A 2 ). Since M evolves the discrete system by one time step, we 
see that 

We iunM ( e ~ik f jAx f + Re ik f jAxA ■ < q 

v ) (A.25) 




is a solution of the finite difference equations V™ + = MVJ 1 across the entire 
grid. Here, the complex frequency uj is given by the dispersion relation (A. 8). 
The solution (A.25) represents a wave that travels from the fine grid region to 
the coarse grid region. At the interface it splits into reflected and transmitted 
pieces. The coefficient R is the reflection coefficient, and T is the transmission 
coefficient. 

The results of this analysis can be checked numerically. For example, in Fig. A.l 
the solid lines show the magnitude of the reflection coefficient calculated from 
Eq. (A. 23). The squares and triangles display the results of a numerical test, 
in which we modeled the propagation of a sine wave contained in a broad 
Gaussian envelope, 

(f)(x, t) = Ae- {x - XQ -^' w2 sin k(x - t) . (A.26) 



The initial data were obtained by discretizing <fi and d<f)/dt at t — 0. The 
constant xq was chosen so that initially the wave packet was situated in the 
fine grid region, far from the interface. We set the Gaussian half-width to 
w = 63Ax, significantly larger than the longest wavelength shown in the 



26 



figures. The magnitude of R was obtained by evolving the wave packet past 
its interaction with the interface, typically about 1000 time steps, and then 
extracting from the numerical data the largest value of \<p\ in the reflected 
pulse. The Courant factor for these numerical runs was At/Axf = 0.4. 

The squares show raw numerical data. The deviation of these data from the 
analytical curves is due to dissipation and dispersion of the wave pulse. Dissi- 
pation can be accounted for rather easily, using the analytical result given in 
Eq. (A. 12). The triangles show the numerical data with correction for dissi- 
pation. The numerical and analytical results for the reflection coefficient are 
in excellent agreement, as seen in Fig. (A.l). 

The above analysis can be repeated for other choices of GCF. In particu- 
lar, for the quadratic GCF of Eqs. (19)-(21), the reflection and transmission 
coefficients are 



Ej(-1 + QEj + 3Ej + £ c 2 (-15 + 10E} - 3Ej)) 
~ 3 + QEj - Ej - £ c 2 (3 - WE} + 15Ej) 

2E c (-3 - 2Ej + IE) + 3E})/E f 
= 3 + QEj -Ej- £ c 2 (3 - lOEj + 15Ej) ' 



with E c and Ef defined as before. Figures (5)-(7) compare the reflection and 
transmission coefficients for linear, direct linear, and quadratic GCF. From 
Fig. 5 we see that quadratic GCF produces a much smaller reflected wave 
than either linear or direct linear GCF. For example, at 20 fine grid points 
per wavelength, linear (direct linear) GCF produces a reflected wave with 
amplitude about 5.6% (4.3%) that of the incident wave. With quadratic GCF 
the reflected wave amplitude is only about 0.31% that of the incident wave. 
For the transmitted wave, Fig. 7 shows that the phase error at 20 fine grid 
points per wavelength is smaller in magnitude by more than a factor of 10 
with quadratic GCF compared to linear or direct linear GCF. Figure 6 shows, 
perhaps surprisingly, that direct linear GCF does the best job of keeping the 
magnitude of T close to 1 at wavelengths less than about 28 Arc/. For longer 
wavelengths, A > 28Axf, quadratic GCF is best. Note, however, that with 
any of these choices of GCF, the deviation of |T| away from unity is relatively 
small. Over the entire range shown in the graphs, A > lOAxj, the maximum 
error for quadratic GCF is less than 1%. For the reflection coefficient, on 
the other hand, the maximum error for direct linear GCF is about 11%. For 
most situations the problem of spurious reflections at an interface will be more 
severe than the problem of inaccurate wave transmission through the interface. 
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A. 2. 2 Discussion of guard cell filling 

The results thus far indicate that quadratic GCF is generally better than 
either linear or direct linear GCF at keeping the reflection coefficient small 
and the transmission coefficient close to unity. The question naturally arises: 
can one do better than quadratic GCF? We will restrict our attention to rules 
for GCF that use a three-point stencil. That is, the fine and coarse grid guard 
cell values V™ and Vq are obtained from linear combinations of three grid 
points, 

Vg = ciVl/2 + c 2 V™ 1/2 + c 3 y" 3/2 , (A.28) 
K = hV?,2 + hV? 1/2 + /3VT3/2 , (A.29) 



where ci, f\, etc. are constants. 

The analysis of Sec. A. 2.1 can be repeated with the guard cells defined as 
above. The resulting reflection and transmission coefficients are functions of 
the constants Ci, f±, etc. Although there are certain combinations of the con- 
stants that outperform quadratic GCF at low resolution, quadratic GCF is 
unique in the following sense. If we consider the high resolution limit, in which 
kfAxf is small, quadratic GCF yields 

R = ^i(k f Ax f f + 0(k f Ax f ) 4 , (A.30) 
T = 1 - ^(k f Ax f ) 3 + Q(k f Ax f r . (A.31) 



The magnitudes of the reflection and transmission coefficients behave like 
|R| = O(kfAxf) 3 and |T| = 1 + O(kfAxf) 41 . For any other choice of con- 
stants Ci, /1, etc., the reflection and transmission coefficients approach and 
1, respectively, more slowly (if at all) than for quadratic GCF. For example, 
for direct linear GCF, we have 

R=U(k f Ax f ) + 0(k f Ax f )\ (A.32) 
T = 1 + ^i(kfAxf) + 0(k f Ax f ) 2 . (A.33) 

Thus, for high resolution, quadratic GCF is the best possible choice given the 
three-point stencil (A.28)- (A.29). Because R and T approach and 1 rapidly, 
as the third power of kfAxf, quadratic GCF performs well at all resolutions. 
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Fig. 1. The evolution of a Gaussian wavepacket on a uniform grid according to the 
linear 1— D wave equation is shown for two different resolutions: H = Ax = 0.045 
(dotted line) and h = H/2 (solid line). 
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Fig. 2. The evolution of the absolute errors e for the two runs in Fig. 1 is shown. 
The dotted line shows en and the solid line 4 x e^. 
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Fig. 3. An interlevel boundary between a coarse and fine grid in 1-D is shown. The 
coarse grid data points are marked with filled diamonds and positive half integers, 
the fine grid data points are marked with filled circles and negative half integers. 
The coarse and fine guard cells are marked with the corresponding open symbols 
and are denoted by "G" and "g," respectively. 
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Fig. 4. The evolution of the absolute errors e is shown for a Gaussian packet crossing 
a fixed refinement boundary at x = 2.1 with linear (dotted line) and quadratic 
(dashed line) GCF. The solid line shows e for a unigrid run at the fine grid resolution. 
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Fig. 5. Absolute values of the reflection coefficient |R| for linear (dashed curve), 
quadratic (solid curve) , and direct linear (dotted curve) GCF for a wave with wave- 
length A crossing a single fixed mesh boundary. The resolution of the fine grid is 
h. 
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Fig. 6. Absolute values of the transmission coefficients for linear (dashed line), 
quadratic (solid line), and direct linear (dotted) GCF. 
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Fig. 7. The phase of the transmission coefficients for linear (dashed line), quadratic 
(solid line), and direct linear (dotted) GCF. 
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Fig. 8. The time evolution of the solution to the nonlinear 1-D wave equation 
is shown. Unigrid runs are shown at the coarse resolution H (thin solid line) and 
the fine resolution h = H/2 (thick solid line). Runs in which the packet encounters 
a mesh refinement boundary at x = 2.1 are shown for linear interpolation (dotted 
line) and quadratic (dashed line) GCF. 
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Fig. 9. Paramesh default linear GCF in 2-D. Points on the coarse grid are denoted 
by diamonds, and points on the fine grid by circles. First, the coarse grid guard cells 
(open diamonds) are filled by averaging the surrounding fine grid points. Then the 
fine grid guard cells (open circles) are filled by taking linear combinations of the 
four surrounding coarse grid points. 
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Fig. 10. Quadratic GCF in 2-D. Coarse grid guard cells (open diamonds) are tem- 
porarily filled by averaging the surrounding fine grid points. Intermediate values 
(open boxes) are obtained by linear interpolation of coarse grid values (solid dia- 
monds) parallel to the interface. Fine grid guard cells (open circles) are then filled 
from a quadratic fit using two fine grid points and one intermediate value. Finally, 
the coarse grid guard cells are filled by matching the coarse and fine grid first 
derivatives across the interface. 
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Fig. 11. The evolution of <fi with a linear 2-D wave equation is shown at 4 consecutive 
times for a run using quadratic GCF to set the data in the guard cells. Each of the 
grid blocks shown has 8x8 zones. 
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Fig. 12. <p is shown along a portion of the x— axis for the 2-D linear wave equation. 
The solid line shows the results of a unigrid run at the fine grid resolution h = 0.0225, 
the dotted line the results of a run with linear GCF, and the dashed line a run with 
quadratic GCF. 
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Linear Interpolation 




Fig. 14. Evolution of gravitational waves in the 3-D Einstein equations using linear 
GCF. g zz — 1 is plotted In the z = plane. Three levels of resolution (h,2h,4h) are 
used, with h = 0.0416667. Each of the grid blocks shown has 6x6x6 zones. 
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Quadratic Interpolation 




15. Same as Fig. 14 except that quadratic GCF is used. 
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Fig. 16. The evolution of g zz — 1 is shown along the y-axis for the 3-D Einstein 
equations. The analytic solution (solid line) and numerical solutions using linear 
(dotted line) and quadratic (dashed line) GCF are shown. The resolution levels are 
given by (h,2h,4h) with h = 0.0416667. 
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Fig. 17. The L2 norm of the absolute error e is shown for two gravitational wave 
runs differing in resolution by a factor of 2 throughout. Both runs have 2 refinement 
boundaries and use quadratic GCF. The solid triangles connected by the solid line 
show e for the run at higher overall resolution, and the solid boxes connected by the 
dotted line show 4 x e for the run with lower overall resolution. 
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Fig. A.l. Absolute value of the reflection coefficient for direct linear guard cell 
filling. The raw numerical data are shown as boxes, while the triangles show the 
data corrected for dissipation. 



46 



